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ABSTRACT 

We examine the evolution and influence of viscosity-induced diskoseismic modes in simulated black 
hole accretion disks. Understanding the origin and behavior of such oscillations will help us to evaluate 
their potential role in producing astronomically observed high-frequency quasi-periodic oscillations in 
accreting black hole binary systems. 

Our simulated disks are geometrically-thin with a constant half-thickness of five percent the radius 
of the innermost stable circular orbit. A pseudo-Newtonian potential reproduces the relevant effects 
of general relativity, and an alpha-model viscosity achieves angular momentum transport and the 
coupling of orthogonal velocity components in an otherwise ideal hydrodynamic numerical treatment. 

We find that our simulated viscous disks characteristically develop and maintain trapped global 
mode oscillations with properties similar to those expected of trapped g-modes and inner p-modes in 
a narrow range of frequencies just below the maximum radial epicyclic frequency. Although the modes 
are driven in the inner portion of the disk, they generate waves that propagate at the trapped mode 
frequencies out to larger disk radii. This finding is contrasted with the results of global magnetohy- 
drodynamic disk simulations, in which such oscillations are not easily identified. Such examples un- 
derscore fundamental physical differences between accretion systems driven by the magneto-rotational 
instability and those for which alpha viscosity serves as a proxy for the physical processes that drive 
accretion, and we explore potential approaches to the search for diskoseismic modes in full magneto- 
hydrodynamic disks. 

Subject headings: accretion, accretion disks-black hole physics-hydrodynamics-X-rays: binaries 



1. INTRODUCTION 

Since the detection of the first high-frequency quasi- 
periodic oscillations ( HFQPOs) froni blac k hole candi- 
date GRS 1915+105 (Morgan et aL]|1997D , much effort 
has been made to relate such oscillations to natural ac- 
cretion disk frequencies. Some early analysis by Nowakj 
et al. ( 1997[ ) suggested that these HFQPOs were man- 
ifestations of global oscillation modes in galactic black 
hole binaries (GBHBs), the theory of which had been ex- 
plored extensively by, for example, Ok azaki et al. (1987) 



and [Nowak fc Wagonerj (|199ll [l992| p^|) I'l 
covery by Strohmayer "([zOUT) of a pair of HFQF 



FQPOs in 

GRO J 165 5-40 with an approximate 3:2 frequency ratio, 
however, lent support to an alternative param etric res- 
onance model (Abramowicz & Kluzniak||2001|) in which 
HFQPOs result from resonance between orbital and ra- 
dial epicyclic motion of disk material. This and similar 
resonance models have the advantage of naturally gen- 
erating the small-integer frequency ratios seen in GRO 
J1655-40 and some subsequently observed sources (for a 
summary of c urrent HFQPO observations, see Remillard 
McClintock|,2006 j . Still, parametric resonance models 
have yet to incorporate convi ncing physical mechanisms 



by which to e xcite HFQPOs (Rebusco 



been noted by |Ortega- Rodriguez et al. (|2008) that mul 



20081), and it has 



tiple global oscillation modes of differing mode number 
produce a 3:2 frequency ratio equally well. A detailed 
physical understanding of such oscillations is crucial since 
their observed frequencies (~ 100 Hz) are comparable to 
orbital frequencies near the innermost stable circular or- 
bit (ISCO) of stellar-mass black holes. Interpreted cor- 



rectly, HFQPOs thus have the potential to tell us much 
about the inner portions of accretion disks and the black 
holes they orbit. 

While numerical simulations have great promise to 
elucidate the nature of disk oscillations, they unfortu- 
nately have failed thus far to produce convincing, iden- 
tifiable HFQPOs. Despite some preliminary claims of 
HFQPO generation in relatively low-resolution simula- 
tions by Kato (2004|), subsequent numerical magnetohy- 
drodynamic (MHDj studies have shown that simulated 
HFQPOs typically are tra nsient ( | Sch nittman et a l.|20Q6l ), 
require external driving (Chan et al. 2006), or, in the 
case of the first paper in this series, remain completely 
undetected ( [Reynolds fc Miller||2008[ hereafter PapefT]). 
In terestingl y, the simulations of ideal hydrodynamic disks 



Paper I did generate trapped gravity- driven {i.e., g- 



mode) global disk oscilla tions such as those described in 



Nowak fc Wagoner ("1992), but these oscillations were not 
seen in their otherwise comparable MHD disks. In fact, 
oscillations of the amplitude seen in X-ray observations 
or in their hydrodynamic disks would have fallen below 
the level of turbulent noise generated in the MHD case, 
so they could not determine whether the modes were hid- 
den or act ively damped by tu rbulence in the manner dis- 



cussed by Arras et al.| (2006). R egardless , it is clear that 
the prototypical MHD disks of [Paper I | failed to excite 
to a detect able level either the globa l diskoseismic oscil- 
lations of Nowak & Wagoner (1993 | or the parametric 
resonance instability of Abr amowicz fc Kluzniak^ (2001). 
In this paper, we describe complementary work to [Pa"^ 



per I in the form of simulations with oscillations inducec 



2 



by the viscous tapping of orbital energy. While the stan- 
dard physical model for black hole accretion is predicated 
upon th e magneto-rotational instability (MRI, "Balbus &' 



Hawley||l991), an intrinsically MHD process that nat- 
urally generates turbulence and the transport of angu- 
lar momentum, the traditional alternative to full MHD 
simulations has been to mimic the influence of magnetic 
flelds t hrough the introduction of an "alpha-model" vis- 
cosity (Shakura & Sunyaev 1973|). This approach sub- 
sumes all physical details of accretion into a single di- 
mensionless viscous parameter a designed to achieve the 
appropriate global level of angular momentum transport. 
This approach is by no mean s completely equivale nt to 
a full MHD t reat ment, as Balbus & Hawley (1998) and 



Pessah et al. (2008) note, and there even remain some ba- 
sic order-of-magnitude discrepancies between values of a 
inferred from observation and those derived from MHD 
simulations |King et al. 2007). Still, studying alpha- 
disks provides us with a method by which to evaluate the 
influence of viscosity independent of MRI-driven turbu- 
lence and other typically complex behaviors of fully MHD 
disks. 



Ortega- Rodriguez & Wagoner (2000) have provided 
a linear analysis of normal modes m viscous, rotating, 
Newtonian fluids that is applicable to our simulated vis- 
cous accretion disks. In particular, they show that the 
presence of viscosity should cause the fundamental g- 
modes in rotating disks to grow at a rate that scales 
with the characteristic orbital frequency in the system. 
In relativistic or pseudo-Newtonian gravitational poten- 
tials, these g-modes are predicted to be non-evanescent 
at radii where |cc;| < z^, where uo is the wave frequency. 



and nis the radial epicyclic frequency [see Nowak fc Wag- 



( 19911 |1992D ; |Ortega-Rodrfguez fc Wagoner] ( |200Uf ; 



or Section 2.2 of [Paper I ' for the appropriate disper- 
sion relations and derivations]. In black hole accretion 
disks, this means that the modes are trapped just under 
the maximum radial epicyclic frequen cy /^max- Similarly, 



Ortega- Rodriguez & Wagoner ( 2000[ ) flnd that viscosity 
also causes the inner pressure-driven (p-mode) oscilla- 
tions to grow for A^max > |^| > ^5 although the successful 
trapping of such modes depends strongly upon on the 

nature of the inner boundary of the disk. 

While the viscosity-induced trapped g-mode of |Ortega^ 
Rodrfeuez & Wagoner] ( |2000D has never been identifled 



explicitly in simulations, numerical models of viscous hy- 
drodynamic disks have generated identiflable waves at 
frequencies comparable to this mode. In an early nu- 
merical analysis of axisymmetr i c, ver tically integrated 
[i.e., ID) disks, 'Honma et al. (1992) found that vis- 
cosity above a critical value a ^ 0.1 caus ed global disk 
oscillations near ^Cmax- Likewise, Chen & Taam (1995) 



and Milsom & Taam ( 1996} identifled global oscillations 
near /^max in vertically integrated disk simulations for a 
range of moderate accretion r ates. This work was fol - 
lowed by the 2D simulations of Milsom & Taam| (1997D, 
which focused on convection in optically thicF7isks but 
also found oscillations near /^max, particularly for lo w ac- 
cretion rate s and large viscosities. More recently, |Mao| 
et al. ([2008) revisited the vertically integrated models of 



Milsom & Taam (1996), pointing out that waves prop- 
agating from the inner portions of the disk could easily 
be locally super-Keplerian. All of these studies associate 
the observed signals with radial inertial-acoustic oscil- 



lations corresponding to the previously mentioned inner 
p- modes. This is a particularly valid interpretation in the 
case of vertically integrated disks where motion is con- 
strained to the radial dimension, but distinguishing be- 
tween trapped g- and inner p-modes in 2D viscous disks 
is not as straightforward, as we shall discuss in this work. 

Our simulations of viscous accretion disks are intended 
to complement this previous work by exploring in more 
detail how viscosity can induce diskoseismic modes in ac- 
cretion disks and how these modes affect the body of the 
disk. First, we seek to discover whether we can produce 
and identi fy in our models any of the v i scosity -induced 



modes of Ortega- Rodriguez fc Wagoner] ( 2000 ) . In par- 
ticular, we are interested in the trapped global g-modes 
since they exist in a narrow frequency range near /^max, 
the value of which in principle can be used as a diagnostic 
of the fundamental physical properties of the black hole. 
Since these g-modes are trapped well away from the inner 
boundary of the disk, we also expect them to be less sus- 
ceptible than inner p-modes to leakage across the ISCO. 
We further examine how such modes generate waves that 
propagate through the entire body of the disk, far beyond 
the formal mode trapping region. Since trapped g-modes 
were identifled in the hydrodynamic simulations - but 
not the MHD disks - of Paper I, we further seek to un- 
derstand and evaluate the observed differences between 
viscous alpha-disks and full MHD models. 

In § [2j we outline the computational framework used 
and describe our simulated disk models. In § [3] we 
present the results of our simulations and discuss our 
identiflcation of trapped diskoseismic modes and the ef- 
fects of these modes on viscous disks. We place our flnd- 
ings in a broader context in § ]4j comparing our results 
to previous work, and present our conclusions in § ]5] 

2. MODELING VISCOUS DISKS 

2.1. Numerical Methods 

To simulate the evolution of viscous accreting systems, 
we have adapted the ZEUS-MP code (version 2), the 



basic workings of which are descr ibed in [Stone fc Nor- 



man 



il992a[li^l 



Stone et al. (1992), and more recently in 
Hayes et al.| ( |2006[ ). This code employs an Eulerian fl- 
nite difference scheme to solve to second-order accuracy 
the equations of ideal compressible fluid dynamics. For 
our purposes, we run ZEUS-MP in pure hydrodynamic 
mode using cylindrical coordinates (r, (j)). The calcula- 
tion is "2.5 dimensional", meaning that it enforces com- 
plete azimuthal symmetry while allowing a non-zero az- 
imuthal velocity. Our simulations feature a gamma-law 
gas equation of state (p ex p^) with a constant 7 = 5/3. 
Zero-gradient outflow boundary conditions are enforced 
at each timestep in both the r and z directions. Addition- 
ally, we employ a protection routine to impose a density 
floor pmin at a value 10~^ times that of the initial disk 
midplane density. 

We have modifled ZEUS-MP to incorporate addi- 
tional physics relevant to the simulation of accretion 
disks. While ZEUS-MP allows the inclusion of point- 
mass gravity, we have adjusted this to reflect a pseudo- 
Newt oniaj]_gravitat^^ such as that devel- 
oped by Paczynski fc Wiita (1980). In this potential. 



GM 



GM 



(1) 
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where R = \Jr^ ^ is the spherical radius. This ap- where r is the cyhndrical radius, z is the vertical height 



proach accurately reproduces the positions of the in- 
nermost stable circular orbit (ISCO) at r = 6rg and 
marginally bound orbit at r = 4rg for a Schwarzschild 
black hole. Additionally, we ha ve added to ZEUS- MP the 
aforementioned "a- model" (Shakura & Sunyaev 1973) 
prescription for viscosity in what is otherwise an ideal 
hydrodynamic system. This modification consists of in- 
troducing a kinematic viscosity of the form 



V — ac.H 



(2) 



where a is a dimensionless constant, Cg is the local sound 
speed, and H ~ c^rjv^ is the scale height of the disk, 
which we can express in terms of the sound speed, cylin- 
drical radius, and azimuthal disk velocity ('U^). This 
model viscosity is applied as a correction to the force 
update in ZEUS-MP and directly updates the velocity 
components exclusively according to 



d{pv) 
dt 



= V • cr. 



(3) 



where the component s of the viscous stress tensor a are 
given in Landau fc Lifshitz ( fl959) , for example. Since 
we include viscosity only as a means to couple velocity 
components and to transport angular momentum, we as- 
sume that the dissipated heat is radiated away instanta- 
neously and thus remove it from the system. For numer- 
ical stability, the viscosity update must take place over 
a timescale less than or equal to the viscous timestep, 
given by 

Advise = C'visc min , (4) 



where Cyisc is a stability constant (Cvisc ~ 0.1) and Axi 
is the computational zone length in the zth dimension. 
In practice, this criterion is met by subcycling the force 
update at a timestep 

Atcour 



At 



sub 



TV 



< At, 



(5) 



where Atcour is the standard Courant-Friedrich s-Levy 
timestep (de scribed in the context of ZEUS-MP by Hayes 
et al. 2006) and N is the smallest positive integer to 
satisfy this condition. 

2.2. Simulated Disk Parameters 

The basic initialization template for our simulated ac- 
cretion disks is identical to that of the 2D hydrodynamic 
disks discussed in (Paper I| Additionally, we conduct two 
"test" simulations (described at the end of this section) 
to confirm that our results do not depend in an unphys- 
ical way upon the details of the initial conditions and 
computational grid size. Since we do not simulate scale- 
dependent processes such as radiative cooling, our disks 
are fundamentally scale-free and we discuss them in nat- 
ural units. 

Our simulated disk density and pressure profiles are 
given by 



p(r, z) = po exp y 



and 



p(r, z) 



GMhl , , 
p(r, z) 



{R-2r^YR 



(6) 
(7) 



above the disk midplane, and R = \Jr^ ^ z^ is again 
the spherical radius. The initial midplane value po is 
independent of radius, as are the scale heights h\ and h^- 
The disk is geometrically thin with a value of /12 = 0.3rg, 
leading to a ratio of /12/nsco = 0.05 at the ISCO. We 
set h\ — I.2/12 so that the disk is ~ 20% too cold to 
maintain vertical hydrostatic equilibrium. As a result of 
this setup, the initial disk collapses and oscillates before 
relaxing into an approximate steady state. The initial 
velocity profile is entirely azimuthal with 



^0 



^fGMT 
r — 2r^ 



Vr = Vz 



0, 



(8) 



for r > nsco- This corresponds to pure Keplerian mo- 
tion in the disk midplane. 

In all simulations, the computational grid spans a ra- 
dial range of r G (4rg,28rg) and a vertical range of 
z e (— 1.5rg, 1.5rg). The grid is populated by zones of 
uniform size Ar ^ 2.3 x 10~^ r^ and Az ^ 1.2 x 10~^ rg, 
leading to a cell aspect ratio of 2:1. This resolution pro- 
vides ~ 25 vertical zones per pressure scale height /12 and 
is thus sufficient to capture waves with wavelengths ~ 4 
times smaller than the scale height. The total duration 
of each simulation is ~ 200Tisco5 where 



Tisco - 61.6 GM/c^ 



(9) 



is the orbital period at the ISCO. 

The only input parameter adjusted across our mod- 
els is the value of the dime nsionless v iscosity param- 
eter a. As summarized in [King et al.| (|2007j), evi- 
dence suggests that observed astrophysical accretmg sys- 
tems feature a ~ 0.1 — 0.4. Rather than restrict our- 
selves to this relatively narrow range of values, how- 
ever, we instead examine an ensemble of simulated 
disks ranging from realistic values to completely in- 
viscid disks. Specifically, we choose model disks with 
a = {0.1, 0.075, 0.05, 0.025, 0.01, 0.0}. This enables 
us to explore how viscosity leads to the development and 
propagation of diskoseismic modes and how this behav- 
ior depends upon the strength of the viscosity. Since the 
value of a is the single criterion by which we distinguish 
our disk models, we refer to them by this value prepended 
with an "A" {e.g., model "A0.05" features a = 0.05). 

Additionally, we briefly describe two "test" models de- 
signed to confirm that the behaviors of our disks are not 
unduly influenced by our specific simulation parameters. 
In the model labeled "EQO.l", we modify the basic disk 
template so that hi = /12 and move the inner edge of 
the disk to a distance r ~ 1.33risco- Since this disk is in 
vertical equilibrium and the inner edge of the disk is com- 
fortably outside of the ISCO, this model helps us to gauge 
how our results depend upon the details of the initial disk 
perturbation. Another model labeled "GRDO.l" features 
a grid that spans a radial range of r G (3.75rg, 28rg) and 
a vertical range of z e (— 1.52rg, 1.52rg). This model 
helps us to identify physically interesting disk behaviors 
and to isolate them from potentially unphysical oscilla- 
tions caused by interactions with the inner radial com- 
putational grid boundary. As their names suggest, both 
test models feature a = 0.1. 

3. SIMULATION RESULTS 
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and f{u) is the Fourier transform 



Fig. 1. — Change of the quantity K = pv^dV with time, 
normahzed to its maximum value. The integration domain V is 
the section r G (7rg,14rg). Shown are models AO (solid), AO.l 
(dotted), AO. 05 (dashed), and AO. 01 (dot-dashed). The long-term 
downward trend in AO shows that the g-modes in this model are 
gradually losing energy that had been provided by the initial disk 
perturbation. The viscous models, on the other hand, rapidly lose 
this initial energy. The evolution of model AO.l, in particular, 
clearly shows that vertical energy is actively replenished in viscous 
disks. 

We now discuss the evolution and analysis of our sim- 
ulated disks. While our models diverge as viscous effects 
become important, each disk is initially perturbed in the 
same way, a nd their e arly behaviors are quite similar. As 
described in [Paper T} the initial disk setup is out of ver- 
tical equilibrium and so falls, rebounds, and eventually 
settles into an approximate steady state. We evaluate 
the decay of the initial disk fluctuations by computing the 
quantity K = pv^dV ^ which is a measure of the energy 
in vertical disk oscillations. Figure [l] shows for models 
AO, AO.l, AO. 05, and AO. 01 the evolution of K in time 
where the integration domain V covers the radial seg- 
ment r G (7rg, 14rg), away from the radial disk bound- 
aries. Although some of the intermediate viscosities are 
omitted from Figure [l] to reduce visual clutter, the mod- 
els shown bracket their behaviors. We will discuss this 
figure in more detail in the following sections, but we 
note here that our different models feature very distinct 
evolutionary profiles in K. This is hardly surprising since 
the viscous runs damp out some of the energy associated 
with the initial disk perturbation while also potentially 
in troducing vertical oscillations through the mechanism 
of Ortega- Rodriguez & Wagoner (2000). Rather than tai- 
lor our analysis of each individual model to that model's 
behavior in we conservatively restrict our discussion 
of all models to times after treiax ~ 6 x 10^ GMjc? ^ cor- 
responding to the approximate exponential decay time 
of K in the inviscid model AO. In [Paper I, this decay 
timescale was seen to increase with higher grid resolu- 
tion, suggesting that numerical dissipation was responsi- 
ble for damping out these initial oscillations. 

3.1. Inviscid Disks 

Let us first review briefly the relevant characteristics 
of AO, the inviscid m odel sim ilar to some of those de- 
scribed extensively in |Paper "T, Since we are most inter- 
ested in physical processes that select specific frequen- 
cies, our primary method of analysis is to compute and 
examine the power spectral density (PSD), defined as 



/(^) = / me 



-27Tiut 



dt 



(10) 



of a given time sequence f{t). Note that the PSDs in 
this paper are taken to be functions of the frequency 
V instead of the angular frequency uo = 27: v. Fig- 
ure [2] shows the midplane {i.e., z = 0) PSDs of the 
radial and vertical velocities, pressure, and density in 
AO. The absolute scales are arbitrary, but one can eas- 
ily see in all four quantities an enhancement approx- 
imately bounded on the right by the radial epicyclic 
frequency. Furthermore, the strongest enhancements 
lie just below the maximum radial epicyclic frequency, 
which, in the Paczynski-Wiita potential, is located at 
^max = cc;Tnax/27r ^ 5. 5 X 10~ ^ c^/GM at a radius of 



7.5 GM/c^. As 



Paper I point out, these features 



have all of the expe cted characteristics of th e trapped g- 



modes described by Nowak & Wagoner '(1992). As noted, 
the evolution of K in Figure [l] illustrates that the energy 
in this trapped mode decays over time, having been intro- 
duced exclusively through the initial disk perturbation. 

Additionally, we note some leakage of the g-mode sig- 
nal both radially inward and outward from rmax- Most 
of the leakage that crosses the IS CO is expected to exit 
the grid, and the leakage radially outward from rmax is 
at least an order of magnitude less in strength than the 
trapped g-mode magnitude for all quantities. We also 
note the presence of some broadband noise near the out- 
ermost disk radii. This signal is caused by disk inter- 
actions with_the outer grid boundary, as the resolution 
tests in Paper I | revealed, and similar features are seen in 
our simulations of viscous disks. 

3.2. Viscous Disks 

Considering now viscous disks, we will focus primarily 
on AO.l for which the a value corresponds most closely 
to dis k viscosities inferred from observations (King et al. 
2007 ). Figure |3] shows the midplane PSDs for the ra- 
dial and verticalvelocities, pressure, and density in AO.l. 
While the colorbars in both Figures [2] and [3] have arbi- 
trary units, they are the same for both figures to fa- 
cilitate cross-comparison of quantities between the two 
models. First, we note that there is extensive broadband 
noise present, particularly for the density and pressure, 
in Figure [3] That this noise is most pronounced at lower 
frequencies is consistent with the prediction by [PapeFT 
that secular variation in the disk caused by the gradual 
loss of material through the radial outflow boundaries 
produces a signal that scales with l/cj^. We have at- 
tempted to remove some of this noise from the density 
and p ressure PSDs using a technique described in |Pa-| 
per I| In this approach, we divide these time series by 
an exponential decay function, choosing the time con- 
stant from a least-square fit to the data. In AO.l, this 
secular trend amounts to a loss of only a few percent of 
the initial total disk mass during the period of analysis 
{i.e., t > treiax), but Figure [s] illustrates that a strong 
residual trend remains. Fortunately, the velocity PSDs 
are affected by this variation only indirectly and, as such, 
require no secular correction. 

The physically significant signal we see in AO.l consists 
of a set of features located near u ^ 5 x 10~^ c^/GM 
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Fig. 2. — Midplane PSDs of radial velocity (top-left), vertical velocity (top-right) , density (bottom-left) and pressure (bottom-right) for 
model AO. Also shown are the radial epicyclic frequency (solid) and orbital frequency (dashed) for comparison. The logarithmic colorbars 
are in arbitrary units and span five or ders of m agnitude. The signal bounded on the right by the radial epicyclic frequency has the properties 
of a trapped g-mode, as described in |Paper~ll 



that extends radially {i.e., vertically in Figure Sl), occu- 
pying up to half of the disk radius. IdentifiaDle in all 
four quantities, this set of features is seen at a frequency 
range very near the maximum radial epicyclic frequency 
of the system at i^max- Although there is some varia- 
tion in the signal profile for the four quantities shown 
in Figure |3j all signals characteristically feature one or 
two prominent broad "spikes" near i^max with multiple 
weaker peaks at adjacent lower frequencies. The spikes 
in the velocity PSDs are seen down to the inner radius of 
the computational grid (r = 4 GM/c^), while those for 
the pressure and density are difficult to identify inward 
of the broadband noise bands near r ~ 7.2 GM/c^. 

To explore the radial dependence of spectral power 
in more detail, we show in Figure [4] a set of PSDs of 
the radial velocity in AO.l linearly added in radius bins 
of width Ar = 1 GM/c^. The vertical bar indicates 
the position of i^max in each plot. At the smallest radii 
(r = 6 — 7 GM/c^), there is sufficient noise in the bin 
to preclude simple identification of the aforementioned 
spikes. While the background noise is slightly dimin- 
ished at r = 7 — 9 GM/c^, the spikes are only easily 



seen by r > 9 GM/c^. Interestingly, the spikes maintain 
a power level that is constant to within a factor of two 
from r = 9 — 15 GM/c^, and they are seen to peak at or 
just below i^max at these radii. In the last three radius 
bins, we see that the amplitude of the signal finally di- 
minishes as lower frequency noise begins to creep in at 
outer disk radii. 

Although PSDs are invaluable for locating such fea- 
tures, we must appeal to alternate methods to determine 
the radius range from which these spikes emanate. Fig- 
urelHshows for model AO.l as a function of both time and 
radius the deviation in midplane radial velocity, defined 
as A^r = Vr — Vr- Showu as a bi-colored dashed line is 
a characteristic outward wave propagation path, derived 
from the local sound speeds averaged over the entire sim- 
ulation time. In this plot, we see that multiple streams 
form a combtooth pattern that, for t > treiax, originates 
at the plateau near r ^ 7.5 GM/c^, runs approximately 
parallel to the dashed line, and finally fades from view by 
r ~ 18 GM/c^. That the pattern runs roughly parallel 
to the dashed line suggests that these are waves moving 
radially outward from their point of origin in the inner 
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Fig. 3. — Midplane PSDs of radial velocity (top-left), vertical velocity (top-right), (decay-corrected) density (bottom-left) and (decay- 
corrected) pressure (bottom-right) for model AO.l. Also shown are the radial epicyclic frequency (solid) and orbital frequency (dashed) 
for comparison. The logarithmic colorbars span five orders of magnitude and are identical to those used in Figure |2] to facilitate cross- 
comparison. The set of vertical spikes just below the maximum radial epicyclic frequency indicate the presence of trapped modes that 
generate waves that pervade the disk. 

disk. Where these waves vanish near the top of Figure 
|5] is near where Figures [S] and [4] suggest that a signal 
associated with the outer grid boundary be gins to m an- 
ifest itself (such features were also noted in Paper I). To 
avoid contamination from such boundary effects, we re- 
strict our analysis to radii inward of r ~ 18 GM/c^, a 
range that still encompasses a large portion of the disk 
external to nsco and rmax- 

We propose that the spikes and associated signals seen 
in Figures [3][5] are the natural result of viscosity-induced 
trapped oscillation mode s, such as those described by 
[Ortega-Rodriguez fc Wagoner (2000). In this framework, 
viscosity provides a mechanism by which the rotational 
velocity of the disk can be tapped by the orthogonal ve- 
locity components, leading to a driven, trapped mode 
(or modes). The narrow frequency range of the spikes, 
located most prominently at or just below i^max, is what 
would be expected from a trapped g-mode, although in 
practice distinguishing between inner p-modes and g- 
modes is not trivial. For example, viscous disks char- 
acteristically share power locally between orthogonal ve- 



locity components, so one cannot simply assume that the 
presence of a signal in indicates a g-mode. Figures [S] 
and|4]in principle could be used to identify the exact ra- 
dial range for these modes, but broadband noise at the 
radii of interest (r ~ 6 — 9 GM/c^) makes it difficult to 
cleanly separate the two distinct mode trapping regions. 
Moreover, it is clear from comparing Figures [3] and [5j for 
example, that one cannot rely upon PSDs to distinguish 
proper trapped modes from induced wave motions. In 
fact, our only strong constraint on the location of these 
modes comes from Figure [Sj which shows that the modes 
themselves are not present exterior to rmax- Given that 
trapped p- and g-modes are apparently indistinguishable 
in our simulations, we will thus refer to these features 
generically as "trapped modes" for the remainder of this 
discussion. 

A combined analysis of Figures |3] - 15] clearly illustrates 
that the influence of these trapped modes in AO.l is not 
restricted to that portion of the disk within the trapped 
region. Figures |3] and [4j for example, demonstrate 
that the disk contains significant power near i^^ax for 
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r > Tmax- Figure [5] shows that this power exterior to rmax 
is in the form of outward propagating waves. That these 
waves are related to the trapped modes is evident from 
their discrete frequencies, which are identical to those of 
the trapped modes. Still, it is challenging to identify ex- 
actly how the modes transfer their energy to these waves 
since both trapped g-modes and inner p-modes are for- 
mally evanescent at frequencies greater than the radial 
epicyclic frequency. Moreover, there is no obvious indica- 
tion that the signal loses power across the radial epicyclic 
boundary, as might be expected for mode leakage. One 



simple plausible explanation is that the modes excite ra- 
dial waves, which are non- vanishing for all f requencies 
greater than the radial epicyclic frequency (see Lubow 



'Pringle| |1993 



for example). One could imagine predom- 
inantly" vertical trapped g-modes, for example, exciting 
through viscous action radial waves that then propagate 
freely outside of the radial epicyclic boundary. Another 
possible explanation is that the trapped modes tunnel 
through the finite evanescent region, which is bounded by 
the radial epicyclic and orbital frequencies in the case of 
the axisymmetric fundamental p-mode. Assuming that 
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Fig. 5. — The deviation in midplane radial velocity (Avr = 
Vt —Vt) for AO.l as a function of radius and time. The linear color 
table extends from Vr = -0.0003c (dark) to Vr = 0.0003c (light), 
where positive radial velocities point radially outward in the disk. 
The bi-colored dashed line, which can be arbitrarily shifted hori- 
zontally, represents the approximate path of a wave moving radi- 
ally through the disk at the local sound speed. That the velocity 
features run parallel to this line suggests that they are waves prop- 
agating radially outward from a region located near t ^max- 

such waves do not decay appreciably in the evanescent 
zone, they could emerge as radial p-modes in regions of 
the outer disk for which the local orbital frequency is 
less than the original trapped mode frequency. Whatever 
the mechanism, the ultimate result is that these waves 
retain the frequency signature of the trapped modes as 
they move through the disk. Eventually, these waves be- 
come lost in the artificial noise generated by the outer 
grid boundary, but not before intersecting a significant 
fraction of the disk body. 

Before moving on to our discussion of other disk mod- 
els, it is worth discussing the slight indications of a signal 
visible near i^max for r < risco in the velocity PSDs in 
Figure [3j This represents leakage of the trapped modes 
down through the IS CO into a narrow accretion stream 
that leaves the grid, similar to that noted in ^Paper I[ 
There is no such leaked signal visible in the density or 
pressure PSDs because these quantities are significantly 
reduced in magnitude radially inward from r^ax 

through 

^isco- The midplane density, for example, is over two or- 
ders of magnitude smaller in the accretion stream than 
at r > Tmax- Any signal proportional to the local den- 
sity would thus have four orders of magnitude less power 
in the PSD at the innermost radii than at r > rmax- 
This typically pushes such a signal below the range of 
our colorbar, and we have confirmed that such signals 
are indeed present, but weak. 

The presence and infiuence of trapped modes are not 
exclusive to model AO.l, but are seen for other viscosi- 
ties as well. Figure |6] shows the radial velocity PSDs and 
Figure[7|the vertical velocity PSDs for AO.l, AO, and the 
four intermediate viscous disk models. The color bars in 
Figures [6] and [7| are each normalized so that we can cor- 
rectly cross-compare magnitudes across different models 
in each figure. First, the basic trapped mode features are 
present in both velocity components for viscous models 
with a > 0.075. While the exact positions and num- 
ber of detectable spikes differ in detail between AO.l and 
AO. 075, the proximity of these spikes to Umax still refiects 
the basic orbital parameters of our model disks. Moving 



to lower viscosities, model AO. 05 shows an identifiable 
signal in vertical velocity that appears to be spatially 
bounded by the radial epicyclic frequency. There is no 
analogous signal in radial velocity, however, even if we 
examine a range in power below that shown in Figure 
[6] Similarly, model AO. 025 shows the hint of a signal in 
vertical velocity for r < rmax, but absolutely no signal 
in the radial velocity. Finally, Model AO. 01 features no 
significant signal for either velocity component or any 
choice of PSD range. 

One might ask whether these trapped modes are iden- 
tical to the trapped g-modes seen in AO or whether these 
are truly the distinct viscosity-induced trapped modes 



described by iOrtega- Rodriguez & Wagoner (2000) 



Looking again at Figure [T] it is clear that the initial disk 
perturbation energy is damped out more rapidly at early 
times in the viscous models than in AO, thus depriving 
alpha-disks of much of the initial energy available to in- 
viscid disks. Moreover, AO shows a long-term downward 
trend upon which is superposed a high-frequency signal 
indicative of the trapped g-mode. Model AO.l, on the 
other hand, has no such obvious long-term trend, sug- 
gesting that viscous and boundary losses are offset by 
ongoing energy input. In the alpha-disk models, this in- 
put energy in fact stems from the viscous coupling of 
disk rotational velocity to radial and vertical motions; a 
channel unavailable to inviscid disks. Although the peaks 
and valleys in Figure T] illustrate that, for model AO.l, 
this process has not achieved a steady-state on timescales 
much shorter than the total simulation time, the overall 
energy profile demonstrates that the viscous method of 
generating persistent trapped modes is distinct from that 
of an inviscid disk. Figure [l] also shows models AO. 05 
and AO.Ol, both of which evolve to a lower value of K. 
This is not surprising since their trapped mode signals 
are weaker or, in the case of AO.Ol, undetectable, sug- 
gesting that the energy resupply is not as efficient as for 
higher viscosities. Interestingly, model AO.Ol does show 
an increase in K near the very end of its evolution, but we 
cannot definitively claim this as evidence for the develop- 
ment of trapped modes without extending the simulation 
in time. 

Finally, we reiterate that the signal observed in our 
viscous models cannot be generated by unphysical com- 
putational phenomena. Test models EQO.l and GRDO.l 
were designed specifically to confirm that conditions such 
as the disk perturbation method and the computational 
grid boundary locations were not important factors in 
our simulations. Figure |8] shows the vertical velocity 
PSDs for the two test simulations. EQO.l resembles 
one of the intermediate viscosity models with a weaker 
trapped mode signal than that of AO.l. This is partly 
because the kinematic viscosity is lower in EQO.l than 
in AO.l. Recall from Equation [2] that u = ac^H ex Cg, 
which is in turn proportional to the disk temperature. 
Since the equilibrium disk in model EQO.l does not col- 
lapse, it is not adiabatically heated and is subsequently 
cooler on average than the disk in model AO.l. Empiri- 
cally, we measure the average kinematic viscosity in the 
inner disk of model EQO.l to be - 70% that of AO.l. We 
thus expect the trapped mode behavior of model EQO.l 
to fall roughly between that of AO. 05 and AO. 075, which 
is consistent with Figure [Sj Additionally, we note that 
the trapped mode signal in EQO.l should take longer 
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Fig. 6. — Midplane PSDs of radial velocity for a range of model viscosities. Also shown are the radial epicyclic frequency (solid) and 
orbital frequency (dashed) for comparison. The logarithmic colorbars span five orders of magnitude and are normalized in magnitude across 
all plots in this figure to facilitate cross-comparison. Here, models AO.l and AO. 075 clearly show trapped modes and waves, and AO features 
a trapped g-mode. 
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studies combine to suggest that oscillations are present 
for a broad variety of physical models, including both 



1996 



Fig. 8. — Midplane PSDs of vertical velocity for the two test simulations EQO.l and GRDO.l. Also shown are the radial epicyclic 
frequency (solid) and orbital frequency (dashed) for comparison. The logarithmic colorbars span five orders of magnitude. EQO.l resembles 
one of the lower viscosity runs because it has a lower average temperature than AO.l. GRDO.l, on the other hand, strongly resembles AO.l. 

than any other model to reach a given trapped mode 
amplitude since the seed disk perturbations are initially 
so small. Model GRDO.l, on the other hand, reflects 
most of the significant characteristics of AO.l, showing 
that the trapped mode signal is not strongly dependent 
in the location of the grid boundaries. Interestingly, this 
and all models feature a small region of high-frequency 
noise near u ~ 0.01 — 0.02 c^/GM. We assume that this 
is associated with motions of material in the relatively 
diffuse accretion stream since the signal is comparable in 
frequency to the innermost orbital frequencies and does 
not appear in the pressure and density PSDs in Figure 



( 19921) , iChen fc Taam 
997p, and |lVlao ef^ 



El 

4. DISCUSSION 

4.1. Comparison with Previous Results 

As mentioned in Secti on [T| t he work of Honma et al. 

19^^ [Milsom Taam| ( |1996[ 
j2008p has identified simulated 
accretion disk oscillations previously, and some of their 
work merits comparison here. Their basic finding rele- 
vant to our work is that waves propagate through their 
simulated disks at frequencies 



constant and alpha viscosities (Milsom & Taam 
1997), and for both 2D and vertically integrated "disEsT 

One interesting feature present in the aforementioned 
viscous disk simulations is a characteristic strong sig- 
nal located in frequency very near i^max- Specifically, 
the higher- frequency spike seen in our models AO.l and 
AO. 075 (see Figures [Sj [4| and [6| peaks at z^max to the 
accuracy of the PSD frequency resolution. This is worth 
noting because analytic treatments of non- evanescent 
trapped g-modes and inner p-modes constrain them to 
have fr equencies strictly less than z/max (e.g., [Orte^a- 
Rodriguez & Wagoner 20001 iNowak & Wagoner || 100:^ 1). 
In fact, the lower- frequency spike seen in models KO.l 
and AO. 075 has exactly the expected characteristics of 
these predicted trapped modes, peaking just below z^max- 
Although the absence of a narrow signal in model AO sug- 
gests that viscosity is partly responsible, no clear physical 
explanation of this high-frequency feature has yet been 
put forth, as 



Kato 



(2001) also notes. 



Specifically, |Mil- 



som & Taam ( 1996 ) find waves in ID vertically integrated 

disks for accretion rates O.OlMsdd ^ ^ ^ 0.25MEdd for 
a characteristic viscous parameter 0.2 < a < 1. For our 
simulated accretion rates of M/M^dd ^ 0.01, we find 
that trapped modes and waves are easily identifiable for 
a = 0.075 — 0.1, trapped modes only are marginally de- 
tectable at a = 0.05, and all modes and waves are un- 
de tectable for lower viscosities. Taking a detailed look 
at Milsom & Taam (1996), we see that they also detect 
a signal for a > 0.05 in the case that M/M^dd = 0.01. 
They claim no detection at a = 0.025, where we too 
failed to detect even a convincing mode, and they simu- 
late no lower viscosities. The 2D simulations of lMilsom fcl 
Taam ( 1997) also find that the oscillations are favored for 
low accretion rates and high viscosities, but it is difficult 
to compare their work to our results since they simulate 
optically thick disks with accretion rates and viscosities 
characteristically higher than ours. We note that these 



4.2. Comparing a-models to Full MHD 

Part of our motivation for conducting these simulations 
was to explore the differences between full MHD and vis- 
cous alpha- models. To d o this, we revisit one of the MHD 
simulations described in Paper I and labeled "MHD_1". 
MHD_1 was a full 3D MHD simulation that utilized a 
computational framework and initial conditions similar 
to our alpha- models, extended axisymmetrically in the 
azimuthal dimension. Additionally, MHD_1 featured ini- 
tially weak poloidal magnetic field loops that threaded 
the disk. As described in [Paper 1\ these fields are am- 
plified by the MRI and drive turbulence which in turn 
provides a natural means for accretion to take place. 
The physical domain of MHD_1 covered r G (4rg, 16rg), 
z G (— 3rg, 3rg),, and (j) G (0,7r/6), and it was run for over 
three times the total simulation time of our alpha-disks. 

Figures [9] and [To] show midplane PSDs of (decay- 
corrected) gas pressure and radial velocity that compare 
MHD_1 (solid), AO.l (dotted), and AO.Ol (dashed). Each 
plot is constructed by summing the power in that quan- 
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Fig. 9. — Midplane PSDs of (decay-corrected) pressure, summed over radial ranges Ar = 0.5rg and each centered at the hsted radius. 
Shown are models MHD_1 (solid), AO.l (dotted), and AO. 01 (dashed). The vertical line indicates the position of the maximum radial 
epicyclic frequency. In all cases, the signals of trapped modes in our model viscous disks would fall at least an order of magnitude below 
the noise level in MHD_1. 



tity over a radial range Ar = 0.5rg centered at the listed 
radius, and the powers are constructed so that cross- 
comparison between models in a given figure is valid. 
In Figure |9j we first notice the remaining secular trend 
in our simulated data that scales in power roughly as 
l/cc;^. The imperfect process of removing this variation 
has left enough of this accretion-related signal present 
that it dominates the overall data trend, particularly in 
the two inner radial bins. On top of this signal, however, 
we do see some oscillations that, in the case of model 
AO.l, are associated with trapped modes. In the three 
bins centered at r > 10 GM/c^, we see two pronounced 
peaks in AO.l that correspond to the waves seen in Fig- 
ures |3" 



[5j for example. At the frequencies of interest near 
lowever, this wave signal is always an order of mag- 
nitude or more below the noise level of model MHD_1, 
which itself shows no convincing trapped mode signal. 
Similarly, there remains no clear trapped mode signal 
in model AO.Ol. Although Figure [lO] features no pro- 
nounced secular trend, we again see that the trapped- 
mode and wave signal in model AO.l is typically more 
than an order of magnitude below the MHD_1 noise level 



at frequencies near i^max- In this case, model AO.Ol does 
show some oscillations near i^max, but these are not obvi- 
ously indicative of trapped mode signals and are always 
at least an order of magnitude below the signal in AO.l. 
Taken together, these two figures suggest that a trapped 
mode signal corresponding to an effecti ve viscosity of the 
magnitude suggested by observations ( King et 31112007 ) 
would not be easily seen above the noise in a real MHD 
disk. 

Our prospects of easily detecting a trapped mode signal 
in MHD_1, however, are even further reduced because of 
that model's low effective viscosity. By measuring the 
correlated stresses, we c an estimate an effective a lpha 
viscosity as described in Balbus & Hawley (1998) and 
|Pessah et al.^ (2008), for example, and given by 



1 



{p5vr5vs) - -—{5Br5Bs) 
47r 



(11) 



Applying this estimator to MHD_1, we find that auuB ^ 
0.01, although the variation in this quantity is compara- 
ble to its value. Still, that auuD <^ 0-1 refiects the 
known discrepancy between simu lated and observa tion- 
ally inferred a values described in King et al. (2007) and 
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Fig. 10. — Midplane PSDs of radial velocity, summed over radial ranges Ar = 0.5rg and each centered at the listed radius. Shown 
are models MHD_1 (solid), AO.l (dotted), and AO. 01 (dashed). The vertical line indicates the position of the maximum radial epicyclic 
frequency. In all cases, the signals of trapped modes in our model viscous disks would fall at least an order of magnitude below the noise 
level in MHD_1. 

further ensures that MHD_1 does not show a peak near 
^max- More reflective of the inferred viscosity in MHD_1 
is the AO. 01 model, for which we have already noted the 
absence of any pronounced trapped mode signal. With 
the addition of turbulent noise at the level present for 
MHD_1, we can safely conclude that this model would 
not generate a detectable signal in a simulated MHD 
disk. 

Although the trapped mode signal present in our 2D 
viscous disks could not be detected in a full 3D MHD 
system, we interpret this more as a shortcoming of the 
approach than as evidence for a dearth of diskoseismic 
modes in real astrophysical syst ems. In addition t o the 
basic discrepancy reported on by King et al. (2007), sev- 



eral groups have begun to address the limit ations of the 
curren t generation of MHD simulations. |Bodo et al.| 
(2008), for example, have noted a dependence upon com- 



putational grid resolution and, surprisingly, grid aspect 
ratio on values of a inferred from shearin g-box sim ula- 
tions of the MRI. Similar l y, bot h Pessah et aLl(|2007) and 
Fromang & Papaloizou (2007) tiave recently discussed 
lack of convergence in zero net magnetic flux shearing- 



box MHD simulations. Specifically, they note that the 
effective alpha viscosities derived from these simulations 
decrease with increasing numerical resolution, suggesting 
that the saturation behavi or of the MRI has yet to be 
captured properly. Pessah et al. I (12007) further point out 
that the physical scales of dissipation in real disks would 
be still smaller than the numerical resolution limit, thus 
making the effective viscosity in analogous real systems 
completely negligible. In these cases, as in the case of 
MHD_1 and the global simulations of Paper I, it is plau- 
sible that field cancellation in the absence of net mag- 
netic flux produces an artificially low effective viscosity. 
Although a strong net vertical field has the p otential to 
disrupt the trapped g-mode region f Fu fc Lai||200 8), this 
process would not affect trapped p- modes. It is also pos- 
sible that one needs only a modest, and therefore non- 
disruptive, net vertical field to seed sufficient turbulence 
to produce a higher effective viscosity. That said, a fac- 
tor of ten increase in a from its inferred value in MHD_1 
would still leave diskoseismic modes at least an order of 
magnitude below the current turbulent MHD noise level, 
making them quite challenging to detect. 
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The presence of turbulent noise in MHD_1 highlights 
one significant way that full MHD simulations are dif- 
ferent from alpha-disks. Our simulated alpha-disks are 
characteristically non-turbulent, particularly as the vis- 
cosity increases. Obviously, this makes the detection of 
diskoseismic modes in alpha-disks simpler because they 
feature less competing background noise than the MHD 
case. This problem can in part be overcome by conduct- 
ing MHD simulations over longer times to produce a bet- 
ter diskoseismic mode signal-to-noise ratio, although one 
must make certain that the integrated mass loss does 
not significantly change the total mass of the disk over 
the simulation time. Such explorations in fact may be 
the only way to correctly ascertain whether the trapped 
modes are hidden beneath the noise or actively damped, 
as suggested by Arras et al. ^ (2006 ) . Finally, we note that 
Pessah et al. ( 2008p have pomted out another shortcom- 
ing of alpha- disks, namely that real MRI-induced stresses 
are not typically proportional to the local shear. All of 
these issues suggest that full MHD treatments are prefer- 
able when net fiux simulations cease to be technically 
prohibitive. 

5. CONCLUSIONS 

We have conducted an ensemble of axisymmetric sim- 
ulations of black hole viscous accretion disks to explore 
the generation of diskoseismic modes and their infiuence 
on disks. While we are still far from a definitive identifi- 
cation of the origin of HFQPOs, we have uncovered and 
explored several interesting facets of viscous disk evolu- 
tion, and we summarize our findings here: 

1) For viscous disks with a > 0.05, we see indications 



of the trapped dis koseismic modes of [Ortega- Rodriguez 



Wagonerl (12000). These modes have all of the expected 



properties of trapped g-modes or inner p-modes, and are 
located at r < rmax with frequencies v ^ z^max- This 



confirms that modes similar to those seen in earlier sim- 
ulations of vertically integrated models are present for 
2D optically- and geometrically-thin accretion disks. 

2) We note that viscous disk models with trapped 
diskoseismic modes also develop related waves that per- 
vade much of the body of the disk. These outward- 
propagating waves are continuous extensions in fre- 
quency and power of the trapped modes, despite extend- 
ing beyond the region of formal mode trapping. This too 
is similar to the ID result and suggests that diskoseismic 
modes can effectively communicate their characteristic 
frequencies to portions of the disk in which the modes 
themselves would be strongly damped. 

3) By comparing our viscous disks to a full 3D MHD 
simulation of Paper I, we have further shown that the 
trapped mode signal for the corresponding alpha disk 
would fall far below the current noise level of the MHD 
simulation. This suggests that, to produce detectable 
trapped modes, MHD simulations may need to feature 
larger effective viscosities, possibly through the natural 
incorporation of net magnetic fiux. Alternately, larger 
trapped mode signal-to-noise ratios should be achievable 
by extending the time domain of these simulations. 
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